Data from the test plate from IMR using “Universal Primers” from Parada et al, 2016. I followed the pipeline from Jesse McNichols and the Fuhrman lab protocols page.

Questions:

Raw Reads

read=read.csv(file='/Users/oliviaahern/Documents/GitHub/TestPlate/raw_seqs.csv',header=T,row.names=1)
data=data.frame(read[,7:8])
#barplot(data$Raw)

No. Reads Community Samples

How many reads in the community samples

sub=subset(read, Frac_Comm=="Comm")
par(mar=c(5,5,0,1))
barplot(sub$Raw,las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,log='x', 
        xlab="Log10 No. Raw Reads")

No. Reads Fractions

sub=subset(read, Frac_Comm=="Frac")
par(mar=c(5,5,0,1))
barplot(sub$Raw,las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,log='x', 
        xlab="Log No. Raw Reads")

No. Reads vs. cDNA Conc.

plot(log10(sub$Raw), sub$Frac.RNA, xlab='Log10 No. of Raw Reads', ylab='cDNA concentration', pch =21, bg='black')

% Proks, Euks, + Cyano

Percent abundance of prokaryotic, eukaryotic, and cyanobacterial reads after sorting raw sequences based on hits to Silva and PR2 databases

On average, 1.8% of the total sequences were Eukaryotes

Community % Proks, Euks, + Cyano

On average 0.46% of the raw reads were eukaryotes, so about 163/47,000 reads.

reads=read.csv(file='/Users/oliviaahern/Documents/GitHub/TestPlate/pct_euk_bact.csv',header=T,row.names=1)
dim(reads)
## [1] 95 17
data1=data.frame(reads[,13:15])

sub=subset(reads, F_C=="Comm")
data1=data.frame(sub[,13:15])


par(mar=c(7,5,1,1),xpd=T)
barplot(as.matrix(t(data1)),las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,
        xlab="% Abundance", col=c('navy','forestgreen','firebrick'),space=0)
box(which='plot')
legend(0,-10, legend=c("Prok", "Cyano", "Euks"), pch =22, pt.bg=c("firebrick", 'forestgreen','navy'),
       bty='n', ncol=3)

Fractions % Proks, Euks + Cyanos

On average, 3.4% of the fractionated sequences were eukaryotes, so about 806/24,000 reads.

Percentage of total raw reads that match to prokaryotic, eukarytoic, or cyanobacterial databases. The samples that don’t reach to 100 did not have matches to the databases, therefore we can assume they are probably chimera.s

sub=subset(reads, F_C=="Frac")
data1=data.frame(sub[,13:15])
barplot(as.matrix(t(data1)),las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,
        xlab="% Abundance", col=c('navy','forestgreen','firebrick'),space=0)
box(which='plot')
legend(0,-10, legend=c("Prok", "Cyano", "Euks"), pch =22, pt.bg=c("firebrick", 'forestgreen','navy'),
       bty='n', ncol=3)

18S sequences

Total of 110 unique ASVs.

read in data

library(phyloseq)
x<-read.csv(file='/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/feature-table.biom.csv',header=TRUE,row.names=1)
OTU = otu_table(x, taxa_are_rows=T)
taxa<-read.csv(file="/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/taxonomy.csv",header=TRUE,row.names=1)
t<-as.matrix(taxa)
tax2<-tax_table(t)
map<-import_qiime_sample_data("/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/map.txt")
#tree=read.tree("/Users/ahern/R/trophic_cascades/mc2/4Nov21/ASVs/tree.nwk")
phyo = phyloseq(OTU, tax2, map)
phyo
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 110 taxa and 84 samples ]
## sample_data() Sample Data:       [ 84 samples by 10 sample variables ]
## tax_table()   Taxonomy Table:    [ 110 taxa by 9 taxonomic ranks ]

overall ASV diversity

library(ggplot2)
#phyo_abund=subset_samples(phyo, Comm_Frac=="Comm")
phyo_abund=transform_sample_counts(phyo, function(x) x/sum(x))
pd <- psmelt(phyo_abund)


colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
#colors=readRDS('colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = ASV)) +
  scale_fill_manual(values=as.character(t(colors))) +
  geom_bar(stat = "identity", width=0.97, color='black',
           lwd=0.1) +
  facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
  labs(x=" ", y = "Class Relative Abundance") +
  theme_bw() +
  theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
        axis.text.x=element_text(size=8, angle=90),
        axis.title = element_text(color='black',face='bold',size=14),
       legend.position = "none")
        
d_rgn

overall ASV diversity w/o Fungi

phy=subset_taxa(phyo, Class!="Fungi")
#phyo_abund=subset_samples(phyo, Comm_Frac=="Comm")
phyo_abund=transform_sample_counts(phy, function(x) x/sum(x))
pd <- psmelt(phyo_abund)




colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
#colors=readRDS('colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
  scale_fill_manual(values=as.character(t(colors))) +
  geom_bar(stat = "identity", width=0.97, color='black',
           lwd=0.1) +
  facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
  labs(x=" ", y = "Class Relative Abundance") +
  theme_bw() +
  theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
        axis.text.x=element_text(size=8, angle=90),
        axis.title = element_text(color='black',face='bold',size=14),
       legend.position = "bottom",
        panel.grid=element_blank(),
        strip.text.x = element_text(
          size = 10, color = "black", face = "bold"),
       legend.text = element_text(size=6),
        legend.key.size = unit(0.25, 'cm')
        #strip.background = element_rect(
        #  color="white", fill="white", size=1, linetype="solid"),
       # panel.spacing = unit(0.05, "lines"))
)
        
d_rgn
## Warning: Removed 1680 rows containing missing values (position_stack).

#sub=subset_samples(phyo_abund, Comm_Frac=="Frac")
#write.csv(otu_table(sub),'18S_nofungi.csv')

#write.csv(tax_table(sub),'18S_nofungi_tax.csv')
sub=subset_samples(phyo_abund, Comm_Frac=="Frac")

control <- subset_samples(sub,SampleID == "Ace.TB5.1" | SampleID == "Ace.TB5.4" | SampleID == "Ace.TB5.5" | SampleID == "Ace.TB5.15" | SampleID == "C12.TB5.7" | SampleID == "C12.TB5.11" | SampleID == "Met.TB5.3" | SampleID == "Met.TB5.9" | SampleID == "Met.TB5.10" | SampleID == "Met.TB5.12")

pd <- psmelt(control)

colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
#colors=readRDS('colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, BD), y = Abundance, fill = Class)) +
  scale_fill_manual(values=as.character(t(colors))) +
  geom_bar(stat = "identity", width=0.97, color='black',
           lwd=0.1) +
  facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
  labs(x=" ", y = "Class Relative Abundance") +
  theme_bw() +
  theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
        axis.text.x=element_text(size=8, angle=90),
        axis.title = element_text(color='black',face='bold',size=14),
       legend.position = "bottom",
        panel.grid=element_blank(),
        strip.text.x = element_text(
          size = 10, color = "black", face = "bold"),
       legend.text = element_text(size=6),
        legend.key.size = unit(0.25, 'cm')
        #strip.background = element_rect(
        #  color="white", fill="white", size=1, linetype="solid"),
       # panel.spacing = unit(0.05, "lines"))
)
        
d_rgn

My T0 compared to Old Siders Data

library(dplyr);library(reshape2);library(ggplot2);library(cowplot);library(tidyverse)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ tibble  3.1.8     ✔ purrr   0.3.4
## ✔ tidyr   1.2.1     ✔ stringr 1.4.1
## ✔ readr   2.1.3     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
theme_set(theme_cowplot())

load("/users/oliviaahern/Downloads/Siders-tagseq-dataframe-26022020.RData",verbose=T)
## Loading objects:
##   siders_tagseq_counts
head(siders_tagseq_counts)
##       DATASET       LOCATION SAMPLEID YEAR REP                       Feature.ID
## 87172   Axial ExtractControl     CTRL 2019     2fa21026b4047e354f119a43c0341a65
## 87236   Axial ExtractControl     CTRL 2019     6bbf315adb7368844fd1dc625be81e30
## 87273   Axial ExtractControl     CTRL 2019     01d46e90107929f3903e4b5bff58a630
## 87616   Axial ExtractControl     CTRL 2019     389950e760558e927385e40108737f9f
## 87957   Axial ExtractControl     CTRL 2019     e613ce213c5added27988e9117e933fd
## 88999   Axial ExtractControl     CTRL 2019     6d4984d3004dc93a767b711fce76c22a
##                                                                                                      Taxon
## 87172                                                                               Eukaryota;Opisthokonta
## 87236                                                                                            Eukaryota
## 87273                                                                                            Eukaryota
## 87616 Eukaryota;Opisthokonta;Fungi;Basidiomycota;Agaricomycotina;Agaricomycetes;Sclerotium;Sclerotium_sp.;
## 87957          Eukaryota;Opisthokonta;Fungi;Basidiomycota;Pucciniomycotina;Cystobasidiomycetes;Rhodotorula
## 88999                                                                                            Eukaryota
##       Confidence    level1       level2 level3        level4           level5
## 87172  0.7041164 Eukaryota Opisthokonta                                      
## 87236  1.0000000 Eukaryota                                                   
## 87273  1.0000000 Eukaryota                                                   
## 87616  0.9742164 Eukaryota Opisthokonta  Fungi Basidiomycota  Agaricomycotina
## 87957  0.9971533 Eukaryota Opisthokonta  Fungi Basidiomycota Pucciniomycotina
## 88999  1.0000000 Eukaryota                                                   
##                    level6                     level7
## 87172                                               
## 87236                                               
## 87273                                               
## 87616      Agaricomycetes Sclerotium;Sclerotium_sp.;
## 87957 Cystobasidiomycetes                Rhodotorula
## 88999                                               
##                               SAMPLE COUNT
## 87172 Axial_ExtractControl_CTRL_2019  1520
## 87236 Axial_ExtractControl_CTRL_2019   260
## 87273 Axial_ExtractControl_CTRL_2019    31
## 87616 Axial_ExtractControl_CTRL_2019   471
## 87957 Axial_ExtractControl_CTRL_2019    24
## 88999 Axial_ExtractControl_CTRL_2019    96
tax=tax_glom(phyo, taxrank="Phyla")
phyo_abund=transform_sample_counts(tax, function(x) x/sum(x))

control <- subset_samples(phyo,Timepoint=="0")

#write.csv(tax_table(control), 'controltax.csv')
#write.csv(otu_table(control), 'controlotu.csv')

pd <- psmelt(control)

#colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
colors <- c("#bfd3e6","#fa9fb5", "#74c476", "#fc8d59", "#807dba", "#238443","#e31a1c", "#ec7014", "#969696")

read=read.csv(file='controlotu.csv',header=T,row.names=1)
colors <- c("#bfd3e6","#fa9fb5", "#74c476", "#fc8d59", "#807dba", "#238443","#e31a1c", "#ec7014", "#969696")

dim(read)
## [1] 9 7
par(mar=c(15,5,1,10),xpd=T)
barplot(as.matrix(read), las=2, col = colors, space =0,
        ylab="% Relative Abundance")
box(which='plot')
legend(8,.7, legend=row.names(read), pch=22,
       pt.bg=colors)

16S Data

Total of 1,632 ASVs.

library(phyloseq)
x<-read.csv(file='/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/16S/feature-table.biom.csv',header=TRUE,row.names=1)
OTU = otu_table(x, taxa_are_rows=T)
taxa<-read.csv(file="/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/16S/taxonomy_16s.csv",header=TRUE,row.names=1)
t<-as.matrix(taxa)
tax2<-tax_table(t)
map<-import_qiime_sample_data("/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/map.txt")
#tree=read.tree("/Users/ahern/R/trophic_cascades/mc2/4Nov21/ASVs/tree.nwk")
phyo = phyloseq(OTU, tax2, map)
phyo
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1632 taxa and 94 samples ]
## sample_data() Sample Data:       [ 94 samples by 10 sample variables ]
## tax_table()   Taxonomy Table:    [ 1632 taxa by 8 taxonomic ranks ]
phyo1 = subset_taxa(phyo, !Order==" Chloroplast")
phyo2 = subset_taxa(phyo1, !Family==" Mitochondria")
phyo=phyo2
phyo2
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1608 taxa and 94 samples ]
## sample_data() Sample Data:       [ 94 samples by 10 sample variables ]
## tax_table()   Taxonomy Table:    [ 1608 taxa by 8 taxonomic ranks ]
agg=tax_glom(phyo2, taxrank="Class")

t=transform_sample_counts(agg, function(x) x/sum(x))
sub=subset_samples(t, Comm_Frac=="Comm")
plot_bar(sub)

phyo_abund=subset_samples(t, Comm_Frac=="Comm")

pd <- psmelt(phyo_abund)


#colors=randomcoloR::randomColor(52)
#write_rds(colors, 'colors_class.rds')
colors=readRDS('/Users/oliviaahern/Documents/GitHub/TestPlate/colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
  scale_fill_manual(values=as.character(t(colors))) +
  geom_bar(stat = "identity", width=0.97, color='black',
           lwd=0.1) +
  facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
  labs(x=" ", y = "Class Relative Abundance") +
  theme_bw() +
  theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
        axis.text.x=element_text(size=8, angle=90),
        axis.title = element_text(color='black',face='bold',size=14),
       legend.position = "bottom",
        panel.grid=element_blank(),
        strip.text.x = element_text(
          size = 10, color = "black", face = "bold"),
       legend.text = element_text(size=6),
        legend.key.size = unit(0.25, 'cm')
        #strip.background = element_rect(
        #  color="white", fill="white", size=1, linetype="solid"),
       # panel.spacing = unit(0.05, "lines"))
)
d_rgn

phyo_abund=subset_samples(t, Comm_Frac=="Frac")

pd <- psmelt(phyo_abund)


#colors=randomcoloR::randomColor(52)
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Frac), y = Abundance, fill = Class)) +
  scale_fill_manual(values=as.character(t(colors))) +
  geom_bar(stat = "identity", width=0.97, color='black',
           lwd=0.1) +
  facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
  labs(x=" ", y = "Class Relative Abundance") +
  theme_bw() +
  theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
        axis.text.x=element_text(size=8, angle=90),
        axis.title = element_text(color='black',face='bold',size=14),
       legend.position = "bottom",
        panel.grid=element_blank(),
        strip.text.x = element_text(
          size = 10, color = "black", face = "bold"),
       legend.text = element_text(size=6),
        legend.key.size = unit(0.25, 'cm')
        #strip.background = element_rect(
        #  color="white", fill="white", size=1, linetype="solid"),
       # panel.spacing = unit(0.05, "lines"))
)
d_rgn
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning: Removed 104 rows containing missing values (position_stack).